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ABSTRACT 



We present 8 years of high-precision radial velocity (RV) data for HD 204313 from 
the 2.7 m Harlan J. Smith Telescope at McDonald Observatory. The star is known 
to have a giant planet (Msini = 3.5Mj) on a ~ 1900-day orbit, and a Neptune-mass 
planet at 0.2 AU. Using our own data in combination with the published CORALIE RVs 
of Segransan et al. (2010), we discover an outer Jovian (Msini = 1.6M/) planet with 
P ~ 2800 days. Our orbital fit suggests the planets are in a 3:2 mean motion resonance, 
which would potentially affect their stability. We perform a detailed stability analysis, 
and verify the planets must be in resonance. 



Introduction 



HD 204313, a sun-like V = 8 star observable in both hemispheres, has been a target for 



multiple RV surveys. Segransan et al. (2010) announced the detection of the first member of 
the star's planetary system with the discovery of HD 204313b, a Jovian-class (Msini ~ 4Mj) 
planet on a long-period (P ~ 5 yr) orbit. More recently, the HARPS survey revealed an interior 
Neptune-mass planet with P = 35 days (Mayor et al. 2011). 



At [Fe/H] = 0.18 (as measured by Segransan et al.|2010[ ), HD 204313 follows the observed trend 



of gas giant hosts being generally metal-rich (e.g. Fischer & Valenti 2005). Furthermore, planet 



c adds to the mounting evidence that Neptune- and lower-mass planets are extremely common 



around main-sequence stars (Howard et al. 2010 Mayor et al. 2011; Wolfgang & Laughlin 2011 



Wittenmyer et al. 2011b). In many ways, HD 204313 represents a "typical" planetary system, 



according to current observations. 



Since 1987, we have used the 2.7 m Harlan J. Smith Telescope at McDonald Observatory for a 
long-baseline RV planet survey (Cochran & Hatzes 1993). An upgrade to our 2D coude spectrograph 



-2 - 



in 1998 gave us access to the full optical wavelength range of our I2 absorption cell, enabling us 
to monitor hundreds of FGK stars with ~ 6 m/s precision over 7-13 years. One of the primary 
scientific objectives of the survey is to obtain a census of Jupiter analogs-giant planets in long period 
orbits (see Wittenmyer et al. 2006, 2011a, for a complete discussion of Jupiter analogs and early 
detection limits from the McDonald Observatory RV survey). We have recently announced three 
giant planets in long-period orbits (Robertson et al. 2012), demonstrating that we have the time 
baseline and sensitivity to detect long-period giants. In the core-accretion theory of giant planet 
formation (Pollack et al. 1996; Lissauer 1995), surface-density enhancement by ices facilitates the 
formation of ~10-15 M ffi cores. The ice line, beyond which ices are present in the protoplanetary 
disk, has been estimated to lie at 1.6-1.8 AU in a minimum-mass solar nebula (Lecar et al. 2006). 



For the case of HD 204313, the inclusion of published CORALIE velocities from Segransan et al. 



(2010) gives us a total time baseline of 12 years, extending our sensitivity comfortably beyond the 
ice line, into the formation locations of gas giant planets. In this paper, we present HD 204313d, 
another Jupiter analog exterior to planet b, and describe its orbital parameters and evolution. 



Observations and Data Reduction 



Our RV data for HD 204313 are all taken from the 2.7 m Smith telescope between July 2003 
and June 2011, resulting in an 8-year time baseline. We use the Tull Coude Spectrograph (Tull et 



al. 1995) with a 1.8" slit, yielding a resolving power R = 60,000. Our RV measurement procedure 



and reduction code AUSTRAL is discussed in detail in Endl et al. (2000). In short, immediately 



before starlight enters the slit, it passes through an I2 absorption cell regulated at 50° C, which 
superimposes thousands of molecular absorption lines over the object spectra in the spectral region 
between 5000 and 6400 A. Using these lines as a wavelength standard, we simultaneously model the 
time-variant instrumental profile and Doppler shift relative to an l2-free template spectrum. The 
resulting RVs are corrected for the motion of the observatory around the solar system barycenter. 
We report our RV data for HD 204313 in Table [1} 



3. Stellar Characterization 



We seek to independently verify the stellar atmosphere parameters for HD 204313 derived by 



Segransan et al. (2010). Using our l2-free stellar template, we measure the equivalent widths of 61 



Fe I lines and 17 Fe II lines. We feed these equivalent widths to the MOOGQ local thermodynamic 
equilibrium (LTE) line analysis and spectral synthesis program ( Sneden||1973| ). By utilizing a grid 
of ATLAS9 model atmospheres ( Kurucz|1993 ), MOOG derives heavy-element abundances to match 
the measured equivalent widths. We then determine effective temperature T e jj by removing any 
trends in abundances versus excitation potential (assuming excitation equilibrium), and computes 
microturbulent velocity £ by eliminating trends with reduced equivalent width (= W\/\). Stellar 
surface gravity is obtained by forcing the abundances measured with Fe I and Fe II lines to match 
(assuming ionization equilibrium). Our measured abundances are differential with respect to the 
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sun. Using a solar port, we have taken a solar spectrum using the same instrumental setup used 
for our RV observations, and run the above analysis for the sun. For reference, we obtain values of 
T eff = 5780 ±70 K, logg = 4.50 ±0.09 dex, f = 1.16 ±0.06 km/s, and loge(Fe) = 7.52 ±0.05 dex 
for the sun. Full details of our stellar analysis can be found in Brugamyer et al. (2011). 

Our stellar parameters for HD 204313 are given in Table [2] For values our routine does not 
calculate, we include catalog values from Kharchenko fe Roeser| ( |2009 Version 3) and Casagrande 



et al. (2011). Our computed values agree extremely well with those presented in Segransan et al. 



(2010). Of particular interest is our measured [Fe/H] of 0.24 ± 0.06, which confirms the metal-rich 
nature of HD 204313. 



4. Orbit Modeling 



Over the 8-year period from July 2003 to June 2011, we have collected 36 RV points for HD 



204313. Our data are plotted as a time series in Figure 1(a) The RMS scatter about the mean 
of these velocities is 40 m/s, with an average internal error of 5.21 m/s. When analyzed alone, 
our RVs show the high-amplitude {K ~ 70 m/s) signal expected as a result of HD 204313b, with 
a period around 5.5 years. We compute Keplerian orbital solutions using the GaussFit modeling 
program ( Jefferys et al.fl9 88) and the SYSTEMIC console ( Meschiari et al.|2009 ), finding excellent 
agreement between the one-planet solutions from both routines. However, a one-planet fit to our 



data gives a period more than 100 days longer than the period reported in Segransan et al. (2010), 



a discrepancy more than 3 times the combined la uncertainties in the orbital period for the two 



models. Additionally, our fit includes a long-period linear trend. We note that Mayor et al. (2011) 



also find a period for planet b considerably longer than the originally published value. 

We then compute a one-planet model using the CORALIE RVs as well as our own. The 
combined RV set includes 132 RVs taken over 10 years. The resulting parameters are closer to the 
previously published solution; we find a period of 2,000 days, with eccentricity 0.16 and a minimum 



mass of 4.36 Jupiter masses. However, this fit is still considerably discrepant from the Segransan 
et al. (2010) solution, and leaves a residual RMS scatter of 11.0 m/s, and a reduced x 2 



5.66. 



We have computed the fully generalized Lomb-Scargle periodogram (Zechmeister Sz Kiirster 



2009) for the combined CORALIE-McDonald data set, and the residual RVs around the one-planet 



fit. The resulting power spectrum is shown in Figure [2} The periodogram of the residual RVs shows 
significant peaks around 340 days, 395 days, and a broad peak between 2700 and 6700 days. We 



calculate a false-alarm probability (FAP) for these peaks using the method described in Sturrock 
&: Scargle (2010), and find a FAP of approximately 5 x 10 -5 for the long-period peak, while the 



395-day and 340-day peaks have FAPs of 4 x 10 4 and 1.5 x 10 3 , respectively. 

We attempt to fit an additional planet to the residuals at each of the periods identified in the 
periodogram. Our fitting routine produces unsatisfactory solutions at the two shorter periods, but 
converges to a fit with an outer giant (Msini = 1.68Mj) planet at 2831 days. The period of planet 
b in this solution is consistent with the original published result, although the eccentricity is higher 
than in a one-planet fit. The resultant two-planet fit is included in Figure l(a)| and we give residual 
plots of the individual planets in Figures 1(b) and 1(c) As with the one-planet solution, GaussFit 



and SYSTEMIC agree nicely on the orbital parameters and their uncertainties. The parameters 
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of our final orbital model are listed in Table |3j The addition of planet d removes the need to 
include a linear slope. Although we attempted to fit an outer planet to each of the CORALIE and 
McDonald RV sets individually, our routines failed to converge for either set. Evidently, both data 
sets are required to achieve the time baseline needed to detect planet d, a fact reinforced by our 
periodogram analysis. When examining the residual RVs to our one-planet fit for each data set 
individually, we see only a monotonic increase in power at long periods for the CORALIE data 
and insignificant power in the McDonald data. Only when the data are combined, and the total 
time baseline exceeds a full orbit of planet d, does the power spectrum show a clearly-defined peak 
around the period of that planet. 



We note that we have not included planet c (Mayor et al. 2011 ) in our analysis. With a reported 



RV amplitude of just 3.28 m/s, the signal of this short-period planet is below the sensitivity limit of 
our data and that of CORALIE. Indeed, our periodogram of the residuals to our two-planet solution 
(Figure]!]) shows no additional signals. Furthermore, the inclusion of a third planet with the orbital 
elements published for planet c does not significantly change our orbital solution. However, the 
RMS (7.79 m/s) and reduced \ 2 (2-98) of our two-planet model are still higher than we expect 
given the precision of the Tull spectrograph and the CORALIE data. While this is reflected in our 
model as a relatively high level of stellar "jitter" (5.46 m/s), our stellar activity analysis (see below) 
suggests HD 204313 should not be so active. Although it would be ideal to include the orbit of 



planet c in our model to verify this hypothesis, we are unable to do so because Mayor et al. (2011 ) 
include neither their measured RVs nor their complete orbital fit for the HD 204313 system. We 
nevertheless conclude that the additional scatter around our fit is most likely due to the unresolved 
planet c, and potentially additional low-mass companions. We refer to the outer planet as HD 
204313d in acknowledgement of the inner Neptune-mass planet. 



Eggenberger et al. (2007) report a companion star 6.2" to the south of HD 204313, although 
they admit a significant probability of a chance alignment. At a distance of 47 pc, the angular 
separation indicates a minimum physical distance of 583 AU between the two objects. The com- 



panion is approximately 9 magnitudes fainter in the near infrared ( Eggenberger et al.||2007 ), and is 



therefore much less massive than HD 204313 if they are in fact bound. If we overestimate the mass 
of this object at 0.5Mq and assume it is associated with HD 204313, the resulting radial velocity 
slope due to the companion is 0.28 m s~ x yr^ 1 , which is roughly equal to our la uncertainty level 
of 0.2 m s _1 yr _1 for a slope in the combined data set. It is therefore safe to conclude that the 
second star is not influencing our modeling of the planetary system. 



5. Stellar Activity and Line Bisector Analysis 

While we do not anticipate that stellar activity should produce RV signals of the amplitudes 
of planets b and d, it is nevertheless important to understand how changes in the atmosphere of 
HD 204313 may influence our velocity measurements, particularly with the amount of scatter seen 
around our fit. We examine stellar activity simultaneously with RV through line bisector analysis 
of stellar lines outside the I2 region and changes in the Ca H and K indices. 

Changes in the stellar photosphere (starspots, etc.) may produce changes in the measured 
RVs. However, these processes will also alter the shapes of the individual stellar absorption lines. 
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Following the method of Brown et al. (2008), we calculate the bisector velocity span (BVS) for each 



of our spectra. The BVS is sensitive to these subtle changes in line shapes, and therefore a reliable 
indicator of activity the stellar photosphere. 

Similarly, if stellar activity is producing RV signals, those signals should also appear in the Ca 
H and K indices. For each RV point in Table [TJ we have computed the Mount Wilson Sjjk index, 



which we list alongside the velocities. From Noyes et al. (1984) we use Shk to derive \ogR' HK , the 



ratio of Ca H and K emission to the bolometric luminosity of the star. From log R'hk we obtain a 
more general idea of the overall activity level of HD 204313. 

All examinations show HD 204313 to be an extremely quiet star. The results of our activity 



analyses are shown in Figure [3j In Figure 3(a) , we plot BVS and Shk versus our measured RVs 
and their residuals around the one-planet fit. In both cases, there is no significant correlation, 
suggesting photospheric activity is not influencing our velocities. Periodograms of Shk and BVS 
(Figure |3(b)[ ) show no periodicity for either index. Furthermore, we measure an RMS of only 17 
m/s for the BVS, and logR' HK = —4.65. It is not surprising, then, that we see no signals or 
correlations in any of our activity indicators. With three planets now known, HD 204313 is rapidly 
becoming a rich planetary system. Its low activity level makes it an ideal candidate for follow-up 
observations to search for additional low-mass companions. 



6. Dynamical Stability Analysis 



A number of recent studies have highlighted the need for observational detections of multiple 
exoplanet systems to be supported by dynamical simulations that test whether the orbits of the 



proposed planets are dynamically feasible (e.g. Horner et al. 2011, 2012a Wittenmyer et al. 2012 



Hinse et al.||2012 ) . Such studies are particularly important when the planets in question appear to 
move on orbits close to mutual mean-motion resonance (e.g. Robertson et al.|2012 ), an architecture 
that can yield either extreme stability or instability, depending on the precise orbits of the planets 
involved. In the case of HD 204313, the best-fit orbits for planets b and d suggest that they may 
well be trapped in mutual 3:2 mean-motion resonance. As such, we chose to perform a highly 
detailed dynamical study of the orbits of planets b and d to investigate whether the orbits that 
best fit the data are dynamically feasible. 



Following earlier work (Horner et al. 2011 Marshall et al. 2010 Robertson et al. 2012), we 



used the Hybrid integrator within the n-body dynamics package MERCURY (Chambers 1999) to 



examine test systems in which the initial orbit of planet b was held fixed at the nominal best fit 
values (in this case, a = 3.04 AU, e = 0.23). The initial orbit of planet d was then systematically 
changed from one simulation to the next, such that scenarios were tested for orbits spanning the 
full ±3(7 error ranges in semi-major axis, eccentricity, longitude of periastron and mean anomaly. 
Such tests have already proven critical in confirming or rejecting planets thought to follow unusual 



orbits (e.g. Horner et al. 2011 Wittenmyer et al. 2012), and allow the construction of detailed 



dynamical maps for the planetary system studied, in orbital element phase space. 

We examined 31 unique values of semi-major axis for planet d, ranging from 3.51 AU to 4.35 
AU, inclusive, in even steps. For each of these 31 initial semi-major axes, we studied 31 values of 
orbital eccentricity, ranging across the full ±3<7 range (e = 0.01 — 0.55). For each of the resulting 
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961 a-e pairs, we considered 11 values of initial longitude of periastron (uj), and 5 values of initial 
mean anomaly (Mo), resulting in a total suite of 52,855 (31 x 31 x 11 x 5) plausible architectures 
for the HD 204313 system. 

In each of these simulations, the masses of the two planets studied were set to their minimum 
(Msini) values. The mass of planet b was therefore set to 3.55Afj, while that of planet d was set 
to 1.68M/. To first order, the more massive the planets, the more strongly they will perturb one 
another, and so setting their masses to the minimum allows us to maximize the potential stability of 
the planetary system. In other words, we expect our resulting dynamical maps to show the maximal 
stability of the orbits tested. The dynamical evolution of the two planets was then followed for a 
period of 100 million years, or until one of the planets either collided with the central star, was 
transferred to an orbit that took it to a distance of at least 10 AU from the central star, or collided 
with the other planet. Collisions were modeled by assuming a density of 1.33 g cm _3 -equal to the 
average density of Jupiter-for each planet and computing a radius accordingly, so that our code 
registered a collision if an actual physical encounter occured. The time of such events was recorded, 



allowing us to construct a dynamical map of the planetary system, shown in Figure 5(a) That 
figure shows the mean lifetime of the HD 204313b-d system as a function of the initial semi-major 
axis, a, and eccentricity, e, of planet d. Each individual initial a-e pair was tested a total of 55 
times, each of which featured a different initial combination of uj-Mq. The lifetimes shown are the 
mean value of the 55 individual lifetimes obtained from those runs. 

Aside from resonant solutions, the entire a-e phase space of allowed orbits for planet d is 
extremely unstable, with collisions between planets b and d often occurring within the first few 
hundred years of the simulations. It is also clear that, even within the resonance, some subset of 
the solutions are dynamically unstable (hence the reason the mean lifetime in the stable region is 
somewhat less than 10 8 years). At the highest eccentricities permitted for the orbit of planet d, no 
stable solutions exist, but there is a broad region of stability within the la errors on the best-fit 



orbit. As can be seen in Figure 5(b) , the stability of orbits in the vicinity of the 3:2 mutual mean- 
motion resonance between the planets is a strong function of the longitude of periastron uj for planet 
d (we note here that planet b's initial longitude of periastron was 298 degrees). Qualitatively, the 
strong uj dependence is reflective of the fact that the 3:2 resonance provides stability by ensuring 
planets b and d never simultaneously approach a true anomaly v ~ 300°, where their orbital paths 
allow very small separations. For configurations outside the stable a — uj space, the resonance 
becomes destructive. Once again, the stable region extends throughout the la uncertainties on the 
best-fit orbit of planet d, with the location of the best-fit orbit lying close to the region of greatest 
stability. 

Here again, our solution suffers from a lack of information regarding planet c. Fortunately, 
at P = 35 days and Msini = 17M©, plausibility arguments suffice to rule out destabilizing 



interactions due to this inner planet. As demonstrated in Horner et al. (2011), planets tend to be 
stable when separated by ~ 5 Hill radii. With planets b and c separated by nearly 12 Hill radii 
(measured from planet b), we expect little mutual influence. A similar argument is presented for 



KOI 961 (Muirhead et al. 2012). Nevertheless, we have performed a small number of simulations 



in which we include a planet with M = 17M®, a = 0.21 AU, and e = 0.17 (Mayor et al. 2011) to 
the stable configurations nearest to our best-fit orbital solution. In all cases, the stability of the 
system is unaffected; while the exact values of a and e for the giant planets are slightly different 
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at each time step when planet c is included, the periods over which the orbital parameters vary 
remain unchanged, and the long-term evolution is the same regardless of whether c is included. We 
therefore conclude that excluding planet c from our larger analysis does not significantly affect our 
results. 

Taken in concert, the results shown in Figure [5] reveal that dynamically stable coplanar solu- 
tions for the orbit of planet d require that it be trapped in mutual mean-motion resonance with 
planet b. Given that the nominal best-fit orbit lies perfectly within the region spanned by that 
resonance, and that a significant fraction of the la error ellipse for planet d is dynamically stable in 
both a-e and a-oj space, we find that our dynamical results are broadly in support of the existence 
of planet d, and may even be used to more tightly constrain its orbit. 



7. Discussion 



With the addition of planet d, HD 204313 joins the growing list of stars hosting multiple gas 
giant planets. At G5 V spectral type and at virtually equal mass to the sun, having two planets 
with masses ~ 2-4 times that of Jupiter makes HD 204313 somewhat of an outlier on the correlation 



between stellar mass and giant planet fraction/mass (Johnson et al. 2011a). However, we confirm 



the star's super-solar metallicity measured by Segransan et al. (2010), thereby offsetting the slight 
discrepancy with stellar mass. 



As of May 2012, there are 12 exoplanet systems in the exoplanets.org (Wright et al. 2011) 
database which contain giant planets believed to be in low-order resonances. However, only HD 



45364 (3:2, |Correia et al.||2009 [) and HD 200964 (4:3, |Johnson et al.||2011b[ ) host multiple gas giants 
in mean-motion resonances closer than 2:1. HD 204313 therefore joins a very small subset of the 
known planet systems. Furthermore, the minimum masses of planets b and d are considerably 
higher than either the HD 45364 or HD 200964 planets, making their continued stability even more 
remarkable. Interestingly, in addition to having the 3:2 resonance in common, HD 45364b/c and 
HD 204313b/d both have mass ratios close to the ~ 3:1 Jupiter/Saturn mass ratio. These systems 



are thus valuable as a comparison to the Nice model (Tsiganis et al. 12005) for the formation of 



the outer solar system. In particular, Batygin et al. (2012) invoke a 3:2 resonance between Jupiter 



and Saturn in simulations which successfully reproduce the orbital configurations of the four outer 
solar system planets and the Kuiper belt. 

The importance of 3:2 mean-motion resonances within the Solar system extends beyond the 
possible interactions between Jupiter and Saturn during the system's early evolution. Beyond the 
orbit of Neptune lie the Plutinos (named after the dwarf planet (134340) Pluto, the first known 
member). These objects, of which several hundred are currently known, are trapped within the 3:2 
Neptunian mean- motion resonance. In Figure [6j we plot the Plutino distribution in a-e space. The 
Plutino population contains objects with a wide range of eccentricities and inclinations, with the 
most eccentric objects crossing the orbit of Neptune, and some moving on orbits that can range as 
close as halfway between the orbits of Uranus and Neptune. The inclinations of the Plutinos range 
from degrees to over thirty degrees. This wide distribution of orbital elements has been used 
to decipher the migration history of Neptune - the idea being that, as that giant planet migrated 
outwards, objects were captured into the 3:2 mean-motion resonance and swept along with the 
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planet, their orbits becoming ever more excited as they were carried along (e.g. Malhotra 1995). 
As a result, the Plutinos nicely map the extent of the stable region of the 3:2 mean- motion resonance 
with Neptune. Though the stable Plutinos do not range to quite as extreme eccentricities as are 
supported for the orbit of HD204313d (as a result of the influence of Uranus on the evolution of the 
most eccentric members), it is striking that the region of stability occupied by the Plutinos is very 
similar to that obtained by our dynamical integrations, as can be seen when comparing Figures 



5(a) and |6| We note in passing that the Hilda family of main belt asteroids are trapped in 2:3 
mean motion resonance with Jupiter, orbiting with periods of ~ 8 years. Despite their sometimes 
high orbital eccentricities (again up to, and in excess of, 0.3), these objects are protected from close 
encounters with the massive planet by the mean-motion resonance they occupy. 



Unlike our recent results for the 2:1 mean-motion resonance in the HD 155358 system (Robert 
son et al.|2012~l ), our dynamical simulations for HD 204313 do not permit coplanar orbits outside the 
3:2 resonance. Evidently, such a small period ratio is only stable when protected by the resonance. 
Furthermore, the existence of unstable orbits within the permitted parameter space shows that 
location within a resonance is not a guarantee of stability. Finally, it is interesting that while for 
both HD 155358 and HU Aquarii ( |Horner et aT|2011[ [Wittenmyer et al.||2012[ ), the inner cutoff for 
stability in a-e space is equal to the inner planet's apastron distance plus five Hill radii, our stable 
solutions for HD 204313 allow planet d to have values of a much smaller than this value. Once 
again, such a conclusion is supported by our knowledge of our Solar system, in which resonant 
configurations ensure the stability of vast populations of objects on orbits that would otherwise 



be highly unstable. Prominent examples include the aforementioned Hildas (e.g. Franklin et al 



1993 


Grav et al. 


2012 


) and Plutinos 


Malhotra 


1995 


Friedland 


2001 


), 


Neptunian Trojans ( 


Morbidelli et al. 


2005 


Sheppard & Trujillc 


2006 



Horner et al. 2012b). In addition to the resonant exoplanets mentioned earlier, these populations 
reinforce the idea that such resonant scenarios are a common outcome from the planet formation 
process. 

We note that the results of our stability analysis not only confirm the validity of our orbital 
model, but in fact place tighter constraints on the system's configuration than our fitting uncertain- 
ties alone. As long-term RV planet surveys such as ours become increasingly sensitive to systems 
with multiple long period companions, it is likely that additional systems with gas giants in close 
resonances will be discovered. Such dynamical simulations are therefore extremely valuable for un- 
derstanding the true architecture of these systems, for which there may otherwise be considerable 
uncertainty as to the orbital parameters. 

It is important to note that our claim that planets b and d are trapped in mutual 3:2 mean- 
motion resonance is based on simulations that assumed that their orbits are coplanar. However, as 
can be seen from the examples of the Hildas and Plutinos within our own Solar system, resonant 
orbits can be dynamically stable for a wide range of mutual inclinations. It might instinctively seem 
that the coplanar case would actually be the least stable configuration, and therefore that mutually 
inclined orbits might allow a broader range of stable solutions. However, we note that in |Horner| 



et al. (2011), the authors considered a wide range of orbital inclinations in an attempt to address 
the apparent instability of the proposed HU Aquarii planetary system, and found that increasing 
the mutual inclination of the planets in question did little to remedy their instability. That said, it 
would certainly be interesting, in future work, to examine the influence of the mutual inclination 
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of the orbits of the planets of the HD204313 system. Fortunately, with predicted astrometric 
displacements of 0.432 mas and 0.264 mas for planets b and d, respectively, both planets should be 
accessible to astrometric measurements with the HST Fine Guidance Sensor (Nelan et al. 2010). 
Plus, the inclusion of inclination constraints would make HD 204313 a unique opportunity for 
comparison to the compact, multi-resonant planet systems discovered by the Kepler spacecraft 
( Holman et al.|2010| pssauer et al.|2011[|Cochran et al.|2011[ ). If astrometry shows the HD 204313 
planets to be coplanar, it would be strongly suggestive that similar migration mechanisms can 
result in systems as different as HD 204313 and the aforementioned Kepler planets. HD 204313 
should therefore be considered a high-priority target for current and future astrometric surveys. 
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Fig. 1. — a. Top: Radial velocity data for HD 204313. Points in black are our 2.7m observations, 



while points in red are CORALIE observations from Segransan et al. (2010). The best-fit orbit 



model is shown as a blue line. Bottom: Residuals to a two-planet fit. b. RVs after subtracting 
our fit to planet d from the velocities in (a). The blue line shows our Keplerian model for planet 
b. c. RVs after subtracting our fit to planet b from the velocities in (a). The blue line shows our 
Keplerian model for planet d. 
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Fig. 2. — From top: [1] Generalized Lomb-Scargle periodogram for the combined 
CORALIE /McDonald RVs of HD 204313. [2] The same periodogram for the residual RVs af- 
ter subtracting a one-planet fit. [3] Periodogram of the residual RVs after subtracting a two-planet 
fit. [4] Periodogram of our time sampling (the window function). The dashed lines indicate the 



approximate power level for a FAP of 0.01, computed from Equation 24 of Zechmeister & Kiirster 
pS09l). 



- 15 - 





(b) 



Fig. 3. — a. Left: Bisector velocity spans plotted against our measured RVs (top) and residual 
RVs to a one-planet fit (bottom) for HD 204313. Right: Shk indices plotted against our measured 
RVs (top) and residual RVs to a one-planet fit (bottom) for HD 204313. 



b. Generalized Lomb-Scargle periodograms for the BVS (top) and Shk indices (bottom) of 
our spectra for HD 204313. The dashed lines indicate the approximate power level for a FAP of 
0.01. 
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Fig. 4. — Face-on orbital diagram of the giant planets in the HD 204313 system. The ellipses shown 
are derived from the model in Table [3] (the open square in Figure [5J , with the lines from the star 
pointing toward the periastron of each planet. The locations of planets b and d are adopted from 
the mean anomalies in Tabled 
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Fig. 5. — a: The mean dynamical lifetime of the HD 204313b-d planetary system, as a function 
of the initial semi-major axis (a) and eccentricity (e) of planet d. The lifetimes are shown on a 
logarithmic scale, ranging from 10 2 years (blue) to 10 8 years (red). The location of the nominal 
best-fit orbit for planet d is denoted by the open square, with the la uncertainties shown by the 
solid lines radiating from that point. It is immediately obvious that the vast majority of the a-e 
space tested is highly unstable, with only a narrow region of stability centered on the mutual 3:2 
mean-motion resonance between planets b and d. 

b: The mean lifetime of the HD 204313b-d system, as a function of the semi-major axis, a, 
and longitude of periastron, u, of planet d's orbit. The lifetime shown at each location in a-u 
space is the mean of 155 individual runs, which tested 31 different orbital eccentricities and 5 
different mean anomalies for that particular a-u) combination. The color scheme is the same as in 
a, and the nominal best-fit orbit for planet d is again marked by the unfilled square, with the la 
errors on that fit shown by the solid lines that radiate from that point. 
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Fig. 6. — The distribution in a-e space of the population of trans-Neptunian objects between 35 
and 43 AU. Note the concentration of objects just beyond 39 AU - the Plutinos, trapped in 3:2 
mean-motion resonance with Neptune. 
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Table 1: Radial Velocities for HD 204313 



BJD - 2450000 


Radial Velocity 


Uncertainty 


Smw 




(m/s) 


(m/s) 




2840.892291 


-8.48 


7.26 


0.1799 ±0.0215 


2933.666733 


-11.81 


4.58 


0.1757 ±0.0231 


3930.951397 


25.41 


4.08 


0.1725 ±0.0262 


4347.740522 


76.44 


6.37 


0.2100 ±0.0327 


4376.744143 


88.50 


4.07 


0.1849 ±0.0217 


4401.617979 


86.26 


4.55 


0.1769 ±0.0250 


4606.945812 


62.97 


5.54 


0.1499 ±0.0233 


4663.924822 


46.99 


5.37 


0.1728 ±0.0251 


4663.938000 


39.25 


2.76 


0.2807 ±0.0466 


4703.749028 


44.40 


5.24 


0.1592 ±0.0206 


4732.798800 


40.00 


3.96 


0.1590 ±0.0214 


4752.702889 


40.71 


6.14 


0.1677 ±0.0217 


4782.617353 


32.02 


6.04 


0.1706 ±0.0228 


4986.953099 


-3.09 


5.05 


0.1604 ±0.0265 


5023.896398 


-10.58 


6.18 


0.1640 ±0.0233 


5049.798314 


-14.30 


6.87 


0.1725 ±0.0246 


5101.700768 


-23.34 


4.92 


0.1682 ±0.0223 


5172.538410 


-20.96 


7.60 


0.1688 ±0.0250 


5470.683390 


-51.54 


4.11 


0.1479 ±0.0222 


5496.568795 


-54.50 


3.56 


0.1583 ±0.0275 


5526.556625 


-54.25 


6.77 


0.2355 ± 0.0439 


5529.556590 


-51.40 


5.55 


0.1704 ±0.0290 


5548.544742 


-46.72 


3.67 


0.1508 ±0.0232 


5722.947515 


-30.38 


7.86 


0.1689 ±0.0294 


5761.927802 


-18.40 


5.41 


0.2000 ± 0.0296 


5790.814330 


-16.89 


4.45 


0.2248 ± 0.0355 


5791.852748 


-16.59 


4.50 


0.1901 ±0.0282 


5811.740128 


-20.36 


6.16 


0.1496 ±0.0251 


5812.750735 


-25.75 


3.56 


0.1498 ±0.0229 


5817.780745 


-33.40 


6.91 


0.2204 ±0.0354 


5838.671924 


-13.10 


5.43 


0.1492 ±0.0231 


5840.686936 


-7.06 


3.68 


0.1768 ±0.0259 


5841.664865 


-10.98 


5.72 


0.1626 ±0.0346 


5842.697923 


-13.76 


5.16 


0.1673 ±0.0367 


5845.589786 


-7.23 


4.97 


0.1803 ±0.0262 


5846.724380 


-18.09 


3.57 


0.2182 ±0.0372 
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Spectral Type 


G5 V 


V 2 


8.006 ± 0.014 


B-V 2 


695 ± 02 


My 


4.63 ± 0.03 


Parallax 2 


21.06 ± 1.04 mas 


Distance 


47 ± 0.3 pc 


T~ f f 


5760 ± 100 K 


log g 


4.45 ±0.12 


[Fe/H] 


0.24 ±0.06 




1.20 ±0.15 km/s 


Mass 3 


1.02M© 


Age 3 


7.20 Gyr 


l°g^W 


-4.65 ± 0.03 



Table 2: Stellar Properties for HD 204313 



Kharchenko & Roeser 



(2009) 



3 From Casagrande et al. (20111, maximum likelihood estimate using Padova isochrones. 



Orbital Parameter 


Planet 


b 


Planet 


d 


Period P (days) 


1920.1 


±25 


2831.6 


± 150 


Periastron Passage T (BJD - 2 450 000) 


2111.6 


±28 


6376.9 


± 176 


RV Amplitude K (m/s) 


57.0 ± 


3 


23.7 ± 


4 


Mean Anomaly Mo 4 


300° ± 


0.4° 


137° ± 


2° 


Eccentricity e 


0.23 ± 


0.04 


0.28 ± 


0.09 


Longitude of Periastron oj 


298° ± 


6° 


247° ± 


16° 


Semimajor Axis a (AU) 


3.04 ± 


0.06 


3.93 ± 


0.14 


Minimum Mass Msini (Mj) 


3.55 ± 


0.2 


1.68 ± 


0.3 


CORALIE RV offset (m/s) 






■19.3 




2.7 m RV offset (m/s) 






29.8 




RMS (m/s) 






7.80 




Stellar "jitter" (m/s) 






5.46 





Table 3: Two-planet orbital solution for the HD 204313 system 



Evaluated at the time of the first RV point reported in Segransan et al. ( 2010 1 



